#ifndef PCPS_PMC_WALKER_HEADER
#define PCPS_PMC_WALKER_HEADER

#include <src/pcps.h>
#include <src/input.h>
#include <src/hamiltonian.h>
#include <src/wavefunction.h>
#include <src/coeff_value.h>
#include <src/move.h>
#include <src/arena.h>
#include <src/random.h>
#include <src/matrix_algebra.h>

namespace pcps { namespace pmc {

  //--------------------------------------------------------------------------------
  // pcps::pmc::Walker -- A class to represent a walker in projector monte carlo
  //--------------------------------------------------------------------------------

  template <class SCALAR, class PART, class REF> class Walker {

    // private data
    private:

      int * _new_config;                            // new lattice configuration
      int * _old_config;                            // current lattice configuration
      double * _matrix_column;                      // normalized column of the transition matrix  I - tau * ( H - E_T )
      double _weight;                               // the weight of the walker
      double _local_energy;                         // local energy for the current configuration
      double _column_norm;                          // normalization constant for the transition matrix column
      int _update_counter;                          // how many times _cvf has been updated without being re-initialized
      pcps::CoeffValueFunc<SCALAR, PART, REF> _cvf; // object for calculating ratios of the trial function's coefficients
      const pcps::action_value_object<SCALAR, pcps::OpRL<SCALAR>,       PART> _avoRL;
      const pcps::action_value_object<SCALAR, pcps::OpNN<SCALAR>,       PART> _avoNN;
      const pcps::action_value_object<SCALAR, pcps::OpN<SCALAR>,        PART> _avoN;
      const pcps::action_value_object<SCALAR, pcps::OpIdentity<SCALAR>, PART> _avoIdentity;

    public:

      // constructor
      Walker(const pcps::Input & userinp,
             const pcps::Hamiltonian<SCALAR> & ham,
             const pcps::Wavefunction<SCALAR> & wfn,
             pcps::Arena & arena)
        : _update_counter(0),
          _cvf(userinp, ham, wfn, arena),
          _avoRL      (userinp.configs_per_site()),
          _avoNN      (userinp.configs_per_site()),
          _avoN       (userinp.configs_per_site()),
          _avoIdentity(userinp.configs_per_site())
      {

        // RLRL operators are not implemented
        if (ham.opsRLRL().size() != 0)
          throw pcps::Exception( (boost::format("pcps::pmc::Walker does not have an implementation for RLRL operators")).str() );

        // allocate space for the lattice configurations
        _new_config = arena.allocate_array<int>(userinp.nsites(), 0);
        _old_config = arena.allocate_array<int>(userinp.nsites(), 0);

        // allocate space for the matrix column
        _matrix_column = arena.allocate_array<double>( 1 + ham.opsRL().size(), 0.0 );

      }

      // functions to access the weight
      double   weight() const { return _weight; }
      double & weight()       { return _weight; }

      // function to access the local energy
      double local_energy() const { return _local_energy; }

      // function to initialize the walker for a given configuration
      template <class GENERATOR>
      void init(const pcps::Input & userinp,
                const pcps::Hamiltonian<SCALAR> & ham,
                pcps::MoveFunc<SCALAR, PART, REF> & mv,
                GENERATOR & gen) {

        // get MPI info
        const MPI::Comm & comm = MPI::COMM_WORLD;
        const int nproc = comm.Get_size();
        const int myrank = comm.Get_rank();

        if (userinp.init_config().size() > 0)
          throw pcps::Exception( (boost::format("pcps::pmc:Walker does not allow user-supplied initial configurations")).str() );

        // set the configuration and initialize coefficient value function object
        _update_counter = 0;
        mv.initialize_configuration(userinp, _cvf, _new_config, gen);
        pcps::bcast(comm, _new_config, userinp.nsites(), MPI::INT, 0);
        for (int i = 0; i < userinp.nsites(); i++)
          _old_config[i] = _new_config[i];
        _cvf.init(_new_config);
        //for (double coeff_norm = 0.0; coeff_norm < 1.0e-12; ) {
        //  mv.initialize_configuration(userinp, _new_config, gen);
        //  try {
        //    _cvf.init(_new_config);
        //  } catch (pcps::SingularMatrixException e) {
        //    coeff_norm = 0.0;
        //    continue;
        //  }
        //  coeff_norm = std::abs(_cvf.initial_value());
        //  for (int i = 0; i < userinp.nsites(); i++)
        //    _old_config[i] = _new_config[i];
        //}

        // compute matrix column and local energy
        this->compute_matrix_column(userinp, ham, true);

        // initialize the weight
        _weight = 1.0;

      }

      // function to copy the state of another walker
      void copy(const pcps::Input & userinp, const pcps::Hamiltonian<SCALAR> & ham, const pcps::pmc::Walker<SCALAR, PART, REF> & other) {
        if (this == &other) return;
        std::memcpy(_new_config, other._new_config, sizeof(int) * size_t(userinp.nsites()) );
        std::memcpy(_old_config, other._old_config, sizeof(int) * size_t(userinp.nsites()) );
        std::memcpy(_matrix_column, other._matrix_column, sizeof(double) * size_t(ham.opsRL().size() + 1) );
        _weight = other._weight;
        _local_energy = other._local_energy;
        _column_norm = other._column_norm;
        _update_counter = other._update_counter;
        _cvf.copy(other._cvf);
      }

      // function to move the walker to a new configuration and update its coefficient accordingly
      template<class GENERATOR>
      void move(const pcps::Input & userinp, const pcps::Hamiltonian<SCALAR> & ham, GENERATOR & gen, const bool fix_nodes) {

        // update the weight
        _weight *= _column_norm;

        // get a random number between zero and one
        double r = pcps::uni_01(gen);

        // choose which column element to move to
        const int col_size = ham.opsRL().size() + 1;
        int i = -1;
        while (r > 1.0e-10 && i < col_size-1)
          r -= std::abs(_matrix_column[++i]);
        assert( i < col_size );

        // make the corresponding move (i == 0 corresponds to staying at the same configuration)
        if (i > 0) {

          // if we are not fixing the nodes, apply the appropriate sign
          if (!fix_nodes)
            _weight *= _matrix_column[i] / std::abs(_matrix_column[i]);

          // get the operator index
          i--;

          // update the new configuration
          _new_config[ham.opsRL()[i].sites()[0]] -= 1; assert( _new_config[ham.opsRL()[i].sites()[0]] >= 0                          );
          _new_config[ham.opsRL()[i].sites()[1]] += 1; assert( _new_config[ham.opsRL()[i].sites()[1]] <  userinp.configs_per_site() );

          // update the coefficient value function object
          _update_counter++;
          //std::cout << boost::format("updating cvf object after moving using operator %i") % i << std::endl;
          if (_update_counter > 1000) {
            _cvf.init(_new_config);
            _update_counter = 0;
          } else {
            _cvf.update(2, &ham.opsRL()[i].sites()[0], _new_config, _old_config);
          }

          // update the old configuration
          _old_config[ham.opsRL()[i].sites()[0]] -= 1; assert( _old_config[ham.opsRL()[i].sites()[0]] >= 0                          );
          _old_config[ham.opsRL()[i].sites()[1]] += 1; assert( _old_config[ham.opsRL()[i].sites()[1]] <  userinp.configs_per_site() );

          // compute the new matrix column and local energy
          this->compute_matrix_column(userinp, ham, fix_nodes);

        }

      }

    private:

      // function to compute the transition matrix column and local energy
      void compute_matrix_column(const pcps::Input & userinp, const pcps::Hamiltonian<SCALAR> & ham, const bool fix_nodes) {

        assert( userinp.pmc_tau() > 0.0 );

        // get the column length
        const int col_size = ham.opsRL().size() + 1;

        // initialize the local energy and matrix column
        _local_energy = 0.0;
        pcps::xscal(col_size, 0.0, _matrix_column, 1);
        _matrix_column[0] = 1.0 + userinp.pmc_tau() * userinp.trial_energy();

        // create an array to hold operators' sites' occupations
        int op_occs[4];

        // add each RL operator's contribution
        for (int i = 0; i < ham.opsRL().size(); i++) {

          // get a reference to the operator
          const pcps::OpRL<SCALAR> & op = ham.opsRL()[i];

          // get the configuration of the ket and check whether it is valid
          bool valid = true;
          for (int k = 0; k < op.sites().size(); k++) {
            int & site_config = _new_config[op.sites()[k]];
            site_config -= op.offsets()[k];
            op_occs[k] = site_config;
            if ( site_config < 0 || site_config >= userinp.configs_per_site() )
              valid = false;
          }

          // if the ket is a valid configuration, process its contribution
          // NOTE:  the sign from fermion operators is taken care of by the get_coeff_ratio function
          if (valid) {

            const SCALAR ratio = _cvf.get_coeff_ratio(2, &op.sites()[0], _new_config, _old_config);

            const double ratio_norm = std::abs(ratio);

            const double x = pcps::real( op.coeff() * _avoRL.get_value(op_occs, userinp.configs_per_site()) * ratio );

            // update the local energy
            _local_energy += x;

            // update the matrix column
            if (x > 0.0 && fix_nodes)
              _matrix_column[0]   -= userinp.pmc_tau() * x;
            else
              _matrix_column[i+1] -= userinp.pmc_tau() * x;
            //_matrix_column[i+1] += userinp.pmc_tau() * std::abs(x);

            //std::cout << boost::format("column element %4i =  %40.12f") % (i+1) % (-userinp.pmc_tau() * x) << std::endl;

          }

          // revert back to the bra configuration
          for (int k = 0; k < op.sites().size(); k++)
            _new_config[op.sites()[k]] += op.offsets()[k];

        }

        // add each NN operator's contribution
        for (int i = 0; i < ham.opsNN().size(); i++) {

          // get a reference to the operator
          const pcps::OpNN<SCALAR> & op = ham.opsNN()[i];

          // get the operator's sites' occupations
          op_occs[0] = _old_config[op.sites()[0]];
          op_occs[1] = _old_config[op.sites()[1]];

          // compute the contribution
          const double x = pcps::real( op.coeff() * _avoNN.get_value(op_occs, userinp.configs_per_site()) );
          _local_energy += x;
          _matrix_column[0] -= userinp.pmc_tau() * x;

        }

        // add each N operator's contribution
        for (int i = 0; i < ham.opsN().size(); i++) {

          // get a reference to the operator
          const pcps::OpN<SCALAR> & op = ham.opsN()[i];

          // get the operator's sites' occupations
          op_occs[0] = _old_config[op.sites()[0]];

          // compute the energy contribution
          const double x = pcps::real( op.coeff() * _avoN.get_value(op_occs, userinp.configs_per_site()) );
          _local_energy += x;
          _matrix_column[0] -= userinp.pmc_tau() * x;

        }

        // add each Identity operator's contribution
        for (int i = 0; i < ham.opsIdentity().size(); i++) {

          // get a reference to the operator
          const pcps::OpIdentity<SCALAR> & op = ham.opsIdentity()[i];

          // compute the energy contribution
          const double x = pcps::real( op.coeff() * _avoIdentity.get_value(op_occs, userinp.configs_per_site()) );
          _local_energy += x;
          _matrix_column[0] -= userinp.pmc_tau() * x;

        }

        // remove numerical zeros from the column
        for (int i = 0; i < col_size; i++)
          if ( std::abs( _matrix_column[i] / userinp.pmc_tau() ) < 1e-10 )
            _matrix_column[i] = 0.0;

        //// if the fixed node on-site potential has caused the probability of staying put
        //// to become negative, set it to zero instead
        //_matrix_column[0] = std::max(_matrix_column[0], 0.0);

        // update the column norm
        _column_norm = 0.0;
        if (fix_nodes)
          for (int i = 0; i < col_size; i++)
            _column_norm += _matrix_column[i];
        else
          for (int i = 0; i < col_size; i++)
            _column_norm += std::abs(_matrix_column[i]);

        // check that the column norm is positive
        if (fix_nodes && _column_norm <= 0.0)
          throw pcps::Exception( (boost::format("projector column norm is not positive with value %.6e") % _column_norm).str() );

        // check that every column element is non-negative
        for (int i = 0; i < col_size; i++)
          if (fix_nodes && _matrix_column[i] < 0.0)
            throw pcps::Exception( ( boost::format("projector column element %i is negative with value %.6e")
                                     % i
                                     % _matrix_column[i]
                                   ).str() );

        // normalize the column
        pcps::xscal(col_size, 1.0 / _column_norm, _matrix_column, 1);

      }

  };

} // end namespace pmc

} // end namespace pcps

#endif
